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Abstract — The nonlinear response of a material to a large elec- 
tric field (steady or pulsed) often determines the ultimate perfor- 
mance of the material for electronics applications. The formalism 
for understanding nonlinear effects in conventional semiconduc- 
tors is well understood. The formalism is less well developed for 
so-called "smart" materials that are tuned to lie close to the metal- 
insulator transition. Here we show that the nonlinear response of 
a strongly correlated electronic material can be calculated with a 
massively parallel algorithm by discretizing a continuous matrix 
operator on the Kadanoff-Baym contour in real time. We bench- 
mark the technique by examining the solutions when the field van- 
ishes and comparing the results to exact results from an equilib- 
rium formalism. We briefly discuss the numerical issues associ- 
ated with the case of a large electric field and present results that 
show how the Bloch oscillations become damped as the scattering 
due to electron correlations increases. 



I. Introduction 

THE problem of the response of materials used in electron- 
ics to large external fields is important from both a theo- 
retical and a practical point of view. On the theoretical side, the 
basic ideas of nonequilibrium statistical mechanics were devel- 
oped over 40 years ago 1 1 1, 1 2 1, but the formalism has not been 
applied to strongly correlated materials except in approximate 
ways. It is interesting to determine exact results for electronic 
systems in an external field which can be used to benchmark 
these approximate techniques. On the practical side, it is often 
the nonlinear behavior of the material or device that determines 
the ultimate performance within electronics. For example, the 
nonlinear current-voltage characteristic of a p-n junction is crit- 
ical for semiconductor-based switching and digital logic, while 
the nonlinear current-voltage characteristic of a nonhysteretic 
Josephson junction allows for digital logic based on rapid single 
flux quantum (RSFQ) ideas |3 1. A "smart" material is a mate- 
rial that can have its properties altered by changing an external 
system variable like pressure, temperature, or a gate voltage. 
The most common devices with tunability are currently based 
on semiconductors or ferroelectrics, but there is increasing in- 
terest in strongly correlated materials near the metal-insulator 
transition, because they might allow for more tunability than 
their semiconducting counterparts. 

The interest in large electric fields arises as the system di- 
mensions shrink onto the nanoscale. When a feature size is on 
the order of 100 nm, a potential difference of 1 V produces an 
electric field of E ^ 10^ V/cm over the feature area. In ad- 
dition, the military is interested in the robustness of devices to 



large pulsed fields that can arise from natural sources like light- 
ning, or from man-made sources like those employed in elec- 
tronic warfare. These high energy-density short-time pulsed 
fields may be difficult to filter out of a device and can cause 
the device to "burn out". 

A "smart material" tuned to lie close to the metal-insulator 
transition is called a strongly correlated material. The name 
arises from the fact that one needs to take into account the 
electron-electron repulsion in determining how the material re- 
sponds to external perturbations. In conventional metals, insu- 
lator, and semiconductors, it is adequate to ignore the mutual 
electron-electron repulsion, and treat all of the electrons as in- 
dependent, moving in an average field created by the other elec- 
trons. This is the regime where band-theory holds. But as the 
electron-electron interactions are made stronger relative to the 
kinetic energy of the electrons, then the electron correlations 
need to be taken into account, implying that one cannot treat 
the other electrons in an averaged way, but one needs to take 
into account where the electrons are and how they move as ev- 
ery other electron moves. This is the regime of strong electron 
correlations, and as the correlations are increased, many mate- 
rials will undergo a metal-insulator transition called the Mott- 
Hubbard transition. 

The simplest model which takes into account strong electron- 
electron correlations is the Falicov-Kimball model |4|. This 
model has two kinds of electrons, itinerant electrons and lo- 
calized electrons. They interact by a Coulomb repulsion when 
they both occupy the same unit cell of the lattice. If the num- 
ber of itinerant electrons plus the number of localized electrons 
is equal to the number of lattice sites, then the system will un- 
dergo a metal-insulator transition as the Coulomb repulsion is 
increased. This model is not appropriate to describe many real 
materials, but its simplicity allows for many exact results to be 
calculated which are vitally important for benchmarking pur- 
poses. 



II. Formalism 

We consider the Falicov-Kimball (FK) model in the presence 
of an external electric field that is spatially uniform, but can 
be time-dependent, and can have an arbitrarily large amplitude. 
The FK model has two kinds of electrons: itinerant electrons 
with creation and annihilation operators c\ and Ci for conduc- 
tion electrons at site i and localized electrons with the corre- 



spending operators and fi. The FK Hamiltonian is 



n 



E 



(1) 

where is the nearest-neighbor hopping matrix, U is the on- 
site repulsion between c and / electrons, /i is the chemical po- 
tential of the conduction electrons and Ef is the site energy for 
the localized electrons. In the simplest case, we ignore the spin 
of the electrons and assume they are spinless. In the calcula- 
tions presented here, we set ^ = U/2, Ef = —U/2, so that 
(etc) = (/t/) = 1/2; this case is called half filling. 

The electric field E(r,t) is described by a vector potential 
A(r, t) in the Landau gauge where the scalar potential van- 
ishes: 

E(M) = -iq^^- (2) 

c ot 

We assume that the vector potential A(r, t) is smooth enough 
in space, that the magnetic field produced by A(r,t) can be 
neglected. 

The electric field is introduced into the Hamiltonian Q by 
the so-called Peierls' substitution |5|, |6| where we neglect in- 
terband transitions because we are considering only a single 
band (the possible dipole transition from a localized electron 
state to a conduction electron state is also neglected; this as- 
sumption may break down if the localized particles are elec- 
trons, but it cannot break down if they are ions as in a binary 
alloy interpretation of the FK model): 



Uj exp 



he 



A(r, t) ■ dr 



(3) 



The Peierls' substitution represents the effect of the line inte- 
gral of the vector potential on the hopping between sites i (at 
position Hi) and j (at position Rj); in this work tij ^ only 
for nearest-neighbor sites i and j. 

For simplicity we shall study the case of a d-dimensional hy- 
percubic lattice in the limit of large spatial dimensions d ^ oo. 
In this limit, the electron self-energy becomes local, which sim- 
plifies both the formalism and the numerical calculations. This 
approximation corresponds to the dynamical mean-field theory 
(DMFT) limit f7]. The simplest electric field is one that lies 
along the unit cell diagonal 1 8 1 : 



A(t)=A(t)(l,l,...,l). 



(4) 



After the Peierls' substitution, the "band-structure" in the elec- 
tric field becomes 
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with a the lattice spacing which we will take to be one. With 
our choice for the electric field along the diagonal, this becomes 
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he 
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with 



(6) 
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and 



Ek 



— ^ > sin k] . 



(8) 



being generalized energy functions and t* is a renormalized 
hopping parameter: t = t* jl^fd in the limit d — > oo |7|; t* 
will be used as our energy unit. 

We find that many quantities we want to determine involve a 
summation over momenta of functions of e and e. These sum- 
mations can be performed more easily by determining a joint 
density of states for the two energies in Eqs. and (|8j; the 
result in the limit of the infinite dimensions |9| becomes: 



1 



■ exp 
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Hence, a summation over an infinite-dimensional Brillouin 
zone can be re-expressed as a two-dimensional Gaussian inte- 
gral. 

In order to solve the many-body problem, we need to de- 
termine the electronic Green's functions in the presence of the 
electric field. The derivation of formulas for these Green's func- 
tions is more complicated than in the absence of a field, because 
there is no time-translation invariance, so the Green's functions 
depend on two different time arguments. Furthermore, since 
the system evolves in the presence of an electric field, there is 
no simple way to relate the quantum-mechanical state at large 
times to the state at small times. Hence, we evolve the system 
forward in time, then we de-evolve it backwards in time, in or- 
der to properly determine its complete time evolution. Since the 
local /-electron number is conserved, and we are not coupling 
the / electrons to the field, the Hamiltonian is a quadratic func- 
tion of the conduction electron operators. This means the time- 
ordered product can be directly evaluated, and relevant func- 
tional derivatives can be taken to determine the Green's func- 
tions. The algebra is somewhat long and will be omitted here. 
The end result is a series of equations for the so-called local 
contour-ordered Green's function, which is defined with two 
time arguments, each one lying on the Kadanoff-Baym con- 
tour (see Fig. □: g''{t,t') = -{i/h){%c^{t)e\{t')), with the 
time-ordering taking place along the contour, the time depen- 
dence of the fermionic operators being determined by the total 
(time-dependent) Hamiltonian in the Heisenberg picture, and 
the angular brackets denoting a weighted trace over all states 
(O) = Tie-'^'^O/Z; the partition function is Z = Tre 



-f3H 



with {3 — 1/T the inverse temperature. In addition, we need to 
define a local self-energy and an effective dynamical 

mean-field A'^(t, t') in analogy with the equihbrium case in zero 
electric field lITUl. itTTl: 

g-{t,t') ^ JdeJ d£p2{£,£) " (9) 

^"(t, t') - .9^„i(t, t') ~ g^-\t, t') - t') (10) 



+ Mg'Llli^-- ^^-u]^\T\t,t') (11) 

= 9tnlit,t') - 9'-\t,t') - X^{t,t'), (12) 



where gQ{e,e,t,t') is the noninteracting Green's function in 
the presence of the electric field |8|, Qimpit, t') is the impurity 
Green's function in zero field, which is equal to gQ{e,e,t,t') 
at £ = £ = 0, and wi is the occupancy of the /-electrons 
[wi — {f^ f)]- The difference of this system of equations from 
the real-frequency case is that these objects are all continuous 
square matrix operators of time (defined on the Kadanoff-Baym 
time contour) rather than being scalar functions of frequency. 
The inverses are all to be interpreted as matrix inverses. 

We will be interested in the so-called lesser Green's function 
g'^ and self-energy I]< in this work. These functions are ex- 
tracted from the contour-ordered objects by fixing the first time 
argument t to lie on the upper real-time piece of the Kadanoff- 
Baym contour and the second time argument t' to lie on the 
lower real-time piece of the Kadanoff-Baym contour 

The system of the equations (I9t-(I12> can be solved by it- 
eration starting from some initial guess for the self-energy 
'S'^{t,t'). From Eq. (|9j, one can find the local Green's func- 
tion g'^~^{t, t'), which allows us to find the effective dynamical 
mean field A'^(t, t') from Eq. (I10> . Then the impurity equation 
[Eq. ( II in allows us to find a new local Green's function, which 
is employed to find a new self-energy from Eq. (I12> . This pro- 
cedure is repeated until the self-energy has converged to a fixed 
point. We call this iterative solution approach the DMFT algo- 
rithm. 
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Fig. 1 

Kadanoff-Baym integration contour for the time variables. 
The time-domain cutoffs are symmetric at itmax. The 

DIRECTION FOR THE INTEGRATION OF THE LINE INTEGRAL IS INDICATED 
BY THE ARROWS. THE DASHED LINE SCHEMATICALLY SHOWS WHERE WE 
TYPICALLY TURN ON THE ELECTRIC FIELD, AS REPRESENTED BY THE 
VECTOR POTENTIAL; IT IS COMMONLY TURNED ON WHEN THE TIME IS 
EQUAL TO ZERO. NOTE THAT FOR THE LESSER FUNCTIONS, WE CHOOSE 
THE FIRST TIME ARGUMENT ON THE UPPER REAL TIME BRANCH, AND THE 
SECOND TIME ARGUMENT ON THE LOWER REAL BRANCH. WHEN THE 
CONTOUR IS DISCRETIZED, WE USE A STEP SPACING OF Ai ALONG THE 
REAL AXIS, AND A STEP SIZE OF 0.05 ALONG THE IMAGINARY AXIS. ALL 
CALCULATIONS PRESENTED HERE HAVE /3 = 1 CORRESPONDING TO 
TWENTY STEPS ALONG THE IMAGINARY AXIS. 



Once the system of equations in Eqs. (I9t-( I12> is solved, then 
we can determine the properties of the system as a function 
of time. It is convenient to describe response functions like 
the self-energy and the Green's function with the relative iioi 
and the average T time variables (the so-called Wigner coordi- 



nates f\T\): 

U,i = t-t', r=^, (13) 

instead of the two times t and t' . In equilibrium, these functions 
are independent of the average time T, and depend only on the 
relative time t^^i- We perform a Fourier transform of the relative 
time to a real frequency, and examine the Green's function and 
self-energy as functions of frequency. In the nonequilibrium 
case, we do a similar thing, performing the Fourier transform 
with respect to the relative time, and examining how the re- 
sponse functions evolve as a function of the average time. For 
example, we have 

S<(c.,r) = J dt,eie'"*-S^ (T+^-f,T-^-fy (14) 

Note that we must have the first time argument t — T + iroi/2 
lying on the upper branch and the second time argument t' = 
T — irei/2 lying on the lower branch of the Kadanoff-Baym 
contour 

Another interesting quantity is the current density that is 
driven by the external electric field: 

JKT) = - ff< e-k, T, T), 

(15) 

with each vector component identical when the electric field 
lies along the diagonal. The magnitude of the total current den- 
sity in the case where the electric field lies along the unit-cell 
diagonal is then 

3{T) = V~dii{T). (16) 

It is well known that the current response to an electric field is 
strange for a perfect conductor that has no electron scattering. 
Indeed, there is an ac response to a dc field due to the lattice 
periodicity, which does not allow the momentum of the elec- 
tron to get too large before an umklapp scattering event occurs 
with the lattice and changes the sign of the momentum. This 
phenomenon is called a Bloch oscillation |13|, |14|, |15|, and 
it should be seen in any material that is free enough of defects 
and other sources of scattering. No conventional metal has ever 
been grown that has small enough scattering to exhibit Bloch 
oscillations. Instead, the scattering occurs so rapidly, that the 
steady-state current is a constant, which increases hnearly with 
the electric field until nonlinear effects take over Bloch oscil- 
lations have been seen in semiconducting heterostructures 1 16l 

Bloch oscillations are also seen in DMFT, with a time- 
independent electric field {E constant and A{T) = -EcT) (HI: 

,(T,..„(£^)/..*<£^,(.), ,,7) 

producing an oscillating current density {p{£) is the noninter- 
acting density of states, which is equal to the integral of p2 over 
£ and f{e) — 1/{1 + exp(/3£) } is the Fermi-Dirac distribution] . 
The frequency of the oscillation is weioch = eE/h and is called 
the Bloch oscillation frequency. We expect these oscillations to 
survive in the presence of scattering if the field is large enough 
that the relaxation time due to scattering is significantly larger 



than the Bloch oscillation period. The frequency of oscillation 
is undoubtedly too high for the Bloch oscillations to be directly 
observed (wBioch ^ 10^^ Hz). 

III. Numerical algorithm 

There are a significant number of numerical issues that need 
to be taken into account to be able to determine the Green's 
functions, and other properties of strongly correlated electrons 
in a large electric field. To begin, the matrix operators are con- 
tinuous operators defined along the Kadanoff-Baym contour, 
and there is no simple way to find their matrix inverse analyti- 
cally. Furthermore, matrix multiplication implies an integration 
over the Kadanoff-Baym contour 

A ■ B{t, t') = ^ dt"A{t, t")B{t", t'), (18) 

which is a complicated line integral in the complex plane (see 
Fig. Our approach to solve this problem is a common nu- 
merical approach — we discretize the Kadanoff-Baym contour 
and evaluate the line integrals as finite Riemann sums over the 
discretized paths. The matrix operators then become finite- 
dimensional square matrices, whose size is equal to the num- 
ber of points used to discretize the Kadanoff-Baym contour 
Once this has been accomplished, then standard LAPACK and 
BLAS routines can be employed to invert and manipulate the 
discretized versions of the matrix operators. One issue that 
needs to be taken into account though is that the inverse of a 
continuous matrix operator satisfies 

/ dt" A-\t,t")A{t" ,t) = 5^{t,t') (19) 

J c 

with 5c the Dirac delta function on the contour The delta func- 
tion is represented by the inverse of the time step used in the 
discretization of the Kadanoff-Baym contour, but one needs to 
note that this time step changes sign on the lower (real) branch 
of the contour, and it becomes imaginary on the vertical piece 
of the contour One needs to properly take this into account 
before using a matrix inversion routine. 

Next, we need to examine the numerical issues arising in 
Eq. (|9j. We evaluate this equation for a given self-energy matrix 
E'^(t, t'). This requires us to choose values of e and e for the 
two-dimensional Gaussian integration, compute the inverse of 
the matrix g§(e, e), subtract the self-energy, and compute a new 
matrix inverse. There is an exact algorithm that allows us to 
directly compute the matrix inverse of g(;(£, e); this arises from 
the equation of motion for the Green's function, from which the 
inverse operator can be directly read off. The only subtlety is 
to ensure that the inverse operator inherits the correct boundary 
condition from the Green's function. This is not so simple to 
carry out, but using techniques like the discretization scheme 
in Negele and Orland |17| and Kamenev 1 18| provides a sys- 
tematic method to directly compute the matrix inverse which 
satisfies the requisite boundary condition and becomes the ex- 
act matrix-operator inverse in the limit where the discretization 
step size goes to zero. The last matrix inversion is a general 
complex matrix inversion, because the self-energy is complex- 
valued, and has no simple symmetries. Hence, it is the most 



"expensive" matrix inversion that needs to be performed. Next 
the matrix elements of the inverse are multiplied by the relevant 
weighting factors for the integration, and finally we accumulate 
the results over all e and e terms that we choose for the two- 
dimensional integral. Since the integral weights are Gaussian, 
it seems reasonable to employ a Gaussian integral scheme for 
choosing the points on the grids and the weights. Unfortunately, 
since each point in the two-dimensional integration requires one 
full matrix inversion, we need to minimize the number of points 
chosen. As a compromise, we use the following scheme: (i) we 
perform the integral using an = 54 Gaussian integration; (ii) 
repeat with an iV = 55 Gaussian integration; and (iii) average 
the results. The iV = 54 case requires 2916 grid points and the 

= 55 case requires 3025 grid points. We choose to average 
these two results, because terms in the Green's function often 
behave like exp(ice), which can be accurately represented by 
the Gaussian integration until c becomes on the order of the in- 
verse of the grid spacing of the Gaussian integration near e = 0. 
Then, the sampling over the discrete points will no longer can- 
cel, and the Gaussian integration will overestimate the value of 
the integral. Since the grid spacing for A^ + 1 points nearly 
interlaces that for A^, the results of the averaging over A^ and 
N + 1 grids produces accurate results for values of c up to two 
times larger than what is possible for either one alone, and in 
the double integral case, it produces a factor of two reduction 
in computation time from using a Gaussian integration scheme 
with twice the number of points. 

This part of the DMFT algorithm, the calculation of the lo- 
cal Green's function from the self-energy, is easily parallelized. 
One simply ships the self-energy matrix, and the energy vari- 
ables e and e to the individual nodes, generates the relevant 
matrix, performs the inversion, and sends the result back to the 
master node for accumulation. Once the local Green's function 
has been calculated, then we proceed with the remainder of the 
DMFT algorithm to determine the new self-energy matrix. This 
part of the code is not so simple to parallelize, because it must 
proceed in a serial fashion. The only possible parallelization 
will occur if we can use SLAPACK routines to distribute the 
calculation of the matrix inverses over a small set of processors. 
To date, we have not included this element in the computation. 

The algorithm given by Eqs. (I9l-( I12> is iterated until it con- 
verges. As the size of the matrices is made larger, by choosing a 
larger maximal cutoff in time for the Kadanoff-Baym contour, 
or by fixing the maximal cutoff time and reducing the step size, 
then the algorithm slows down significantly, and it becomes 
more difficult to attain the same level of accuracy at the end 
of the iterations. We usually try to iterate the solutions at least 
20 — 50 times for full convergence, but sometimes we have to 
limit ourselves to about 10 iterations due to the computational 
time involved. 

There are a number of numerical issues that play a role in the 
quantitative accuracy of the results. These include the maximal 
time chosen for the cutoff imax, the step size in real time At, 
and the number of points A^ chosen for the e grids. In this work, 
we examine a number of different choices for these parameters 
to see which terms are the most important for maintaining ac- 
curacy of the results. We do this employing what is probably 
one of the most difficult cases for the nonequilibrium code, that 
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Fig. 2 

Lesser SELF-ENERGY FOR /3 = 1, (7 = 1, and half-filling. We 

EXAMINE THE FREQUENCY DEPENDENCE OF THE SELF-ENERGY AT AN 
AVERAGE TIME T = AND IN EQUILIBRIUM, BUT CALCULATED WITH THE 
NONEQUILIBRIUM FORMALISM. THE EXACT RESULTS ARE IN BLUE, AND 
THE OTHER CURVES ARE ALL CALCULATED WITH At = 0.05 AND 

fjnax = 15. The BLACK CURVE USES Gaussian INTEGRATION WITH 

AT = 54 AND 55 POINTS. THE RED CURVE IS FOR = 100 AND 101 
POINTS, AND THE GREEN CURVE IS A TRAPEZOIDAL RULE WITH 1000 
POINTS EVENLY SPACED BETWEEN —3 AND 3. 
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Fig. 3 

Imaginary part of the lesser self-energy for fi = \,V = \, and 

HALF-FILLING. WE EXAMINE THE RELATIVE TIME DEPENDENCE OF THE 
self-energy AT AN AVERAGE TIME T = AND IN EQUILIBRIUM, BUT 

CALCULATED WITH THE NONEQUILIBRIUM FORMALISM. THE EXACT 
RESULTS ARE IN BLUE, AND THE OTHER CURVES ARE ALL CALCULATED 
WITH Ai = 0.05 AND tmax = 15. THE COLORS OF THE OTHER CURVES 
ARE IDENTICAL TO THOSE IN FiG. 2. iNSET IS A BLOW UP OF THE REGION 
AROUND t^gi = -20, WHERE THE INTEGRATION OVER £ INTRODUCES 
SPURIOUS RESULTS WHEN THE GRID SPACING IS TOO COARSE, WHICH 
LEAD TO THE OSCILLATIONS SEEN IN FiG. 2. 



of a vanishing electric field (corresponding to an equilibrium 
situation). There is a major simplification of this approach be- 
cause we have no e dependence in our formulas, and the inte- 
gral over e can be done trivially. The problem arises from the 
fact that the self-energy develops a delta function in frequency 
at w — in the insulating state, and this function cannot be 
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Fig. 4 

Lesser self-energy for /3 = i, !7 = i, and half-filling for 

DIFFERENT DISCRETIZATIONS OF THE KADANOFF-BAYM CONTOUR. THE 
EXACT RESULTS ARE IN BLUE, AND THE OTHER CURVES ARE ALL 
CALCULATED WITH Af = 54 AND 55 POINTS FOR THE GAUSSIAN 
INTEGRATION AND tmax = 15. THE BLACK CURVE HAS At = 0.1. THE 
RED CURVE HAS At = 0.075 AND THE GREEN CURVE HAS At = 0.05. 
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Fig. 5 

LESSER Green's function for /3 = i, f7 = i, and half-filling for 

DIFFERENT DISCRETIZATIONS OF THE KADANOFF-BAYM CONTOUR. WE 
PLOT THE RESULTS FOR T = AS A FUNCTION OF tjoi . THE MAIN PANEL 
IS THE IMAGINARY PART, AND THE INSET IS THE REAL PART. THE EXACT 
RESULTS ARE IN BLUE, AND THE OTHER CURVES ARE ALL CALCULATED 
WITH Af = 54 AND 55 POINTS FOR THE GAUSSIAN INTEGRATION AND 
tmax = 15. THE BLACK CURVE HAS At = 0.1. THE RED CURVE HAS 
At = 0.075 AND THE GREEN CURVE HAS At = 0.05. NOTE THAT THE 
EXACT RESULT IS AN EVEN FUNCTION OF t^cl FOR THE IMAGINARY PART 
AND AN ODD FUNCTION FOR THE REAL PART. 



easily represented in a calculation for real time that has a finite 
cutoff along the time axis. Indeed, we find that it is difficult 
to obtain good results for the self-energy as a function of fre- 
quency in this case. But we also check the moments of the 
self-energy and the Green's function, and find good agreement 
for the low-power moments, indicating that the results do work 
well for short times. Hence, it is reasonable to think that they 
do a good job at determining the initial transient response to an 
electric field after the field is turned on. Furthermore, we do 
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Fig. 6 

Current density as a function of time for U = 0.5 and various 

VALUES of At (tmax = 15 AND = 54, 55 GAUSSIAN INTEGRATION). 

The temperature is fixed at /3 = 1 and the electric field 
SATISFIES E = I. Note how the current converges to a "fixed 

POINT" as At — > 0, BUT NOT UNIFORMLY. 
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Fig. 7 

Current density as a function of time for U = 0, 0.5, and 1. The 
temperature is fixed at 13 = I and the electric field satisfies 
E = 1 (N = 5A, 55 Gaussian integration). Note how the Bloch 
oscillations are damped as the scattering increases, and how 

THE period might BE MODIFIED. 



not know whether the delta function in the self-energy survives 
when a field is turned on. If it doesn't, then the nonequilibrium 
approach may work well for the insulating cases, even if it has 
problems for the equiUbrium system. 

IV. Benchmarking the equilibrium results 

The electronic density of states is a quantum-mechanical 
measure of the elementary excitations of a strongly correlated 
material. In most materials, the density of states spans a 
range of a few electron volts, or in our units, a few t* . The 
Fourier transform of the density of states is closely related 
to the Green's functions in the time domain. We expect the 



time-dependent Green's functions to have important time de- 
pendence when iici is smaller than c/t* for some constant c 
of order one. Hence choosing a cutoff in time on the order of 
10 — 20 1/t* is a reasonable choice. The above heuristic dis- 
cussion holds when the system is metallic. When it becomes 
insulating, then the density of states has a gap, and the self- 
energy develops a pole, for small frequencies. In this case, we 
need a large time domain (infinite in the case of a delta func- 
tion) to properly reproduce the density of states or self-energy 
after Fourier transforming. This means that the time-domain 
cutoff needs to be made larger as the electron correlations are 
made stronger. 

Similarly, the more fine-structure that is present in the 
Green's functions as functions of time, the smaller the step 
size At in the time domain needs to be to be able to accurately 
discretize the matrix operators. In general, we expect the dis- 
cretization error to grow as the electron correlations increase, 
so At will need to be reduced for larger values of U. We typi- 
cally choose At = 0.1 ~ 0.05 in our calculations. Obviously, 
we cannot keep increasing the time domain, and decreasing the 
step size when we have finite computational resources. The 
maximal matrix that we consider has a size on the order of 
2500 X 2500. This choice is not a function of memory limita- 
tions, but rather is an issue of the computational time needed to 
invert these matrices as part of the DMFT algorithm. The bot- 
tom line is that these calculations will not be able to be pushed 
to too large a value for U without running into resource prob- 
lems. 

In order to benchmark our code, we have chosen to examine 
the equilibrium solutions using the nonequilibrium formalism. 
This is a particularly nice exercise to undertake, because the 
equilibrium solutions are all known to high accuracy via a di- 
rect solution using an algorithm in the frequency domain 1 10|, 
ml . It is also a challenging test of the nonequilibrium codes, 
because we need to Fourier transform the solutions in time to 
functions of frequency, and the effects of the discretization step 
size, and of the time-domain cutoff can play critical roles in the 
accuracy. The other important parameter in governing the ac- 
curacy is the step size employed in the energy integrations to 
determine the local Green's function. As discussed above, if 
the step size is too large, then we can generate spurious signal 
at large relative times, which will Fourier transform into high- 
frequency wiggles in the frequency domain. If we know that 
such structures are not present in the exact results, we can eas- 
ily filter them out, but this becomes problematic when we are 
not sure whether such structures are real or numerical artifacts 
(which occurs when we perform nonequilibrium calculations). 

We benchmark our results by examining the equilibrium so- 
lutions at high temperature (/3 = 1) and for a large value of U 
(U ~ 1) that is still in the metallic phase (the metal-insulator 
transition occurs at U = \/2). We choose this for our initial 
benchmarking exercise because the self-energy does not have a 
pole; hence, the numerical issues should be under better con- 
trol. We will briefly discuss issues that occur in the insulating 
phase below. 

Our first result is shown in Fig. |2] It plots the lesser self- 
energy at T = as a function of ui, which is calculated by per- 
forming the Fourier transformation with respect to trci- This re- 



suit should be independent of T, because we are in equilibrium, 
but the results do have a small dependence on T, that is due to 
the fact that we have instituted a finite cutoff in time tmax = 15. 
In this figure, we study how sensitive the results are to changing 
the number of points in the integration over e (recall, the inte- 
gration over e is trivial when we are in equilibrium). The blue 
curve is the exact result, the black curve employs Gaussian in- 
tegration, averaging over the = 54 and iV = 55 cases, while 
the red curve is similar with = 100 and = 101 points. 
The green curve uses a much smaller step size in e, employ- 
ing — 1000 points in a trapezoidal rule integration, ranging 
from —3 to 3. Note how all of the results lie on top of each other 
for small lo, although they do differ from the exact result. In the 
tails, for larger |cj|, we see oscillations develop for the Gaussian 
integrations, that are reduced as the step size is made smaller, 
and completely disappear by the time N = 1000. These results 
show that by carefully controlling the step size used for the en- 
ergy integration, one can obtain converged results without any 
extraneous oscillations, but those converged results are not ex- 
act, because they were calculated with finite values for At and 

In Fig. |3j we show the results for the imaginary part of the 
lesser self-energy as a function of t^^i- Once can see partic- 
ularly good quantitative agreement for small times, and then 
there is a region with oscillations out in the tail (trci ~ —20). 
This region is blown up in the inset. The oscillations are only 
present for the coarse e integrations, and the amplitude of the 
oscillations shrinks as the step size is made smaller. It is pre- 
cisely these spurious results that lead to the oscillations in the 
Fourier transform (see Fig. If we know about this kind of 
spurious behavior, we can filter the oscillations out before per- 
forming the Fourier transform, but this is an ad hoc procedure 
that cannot be generalized to the nonequilibrium case. 

In Fig. |4] we show the At dependence of the calculations 
with fixed values for tmax = 15, and A^ = 54 and A^ = 55 av- 
eraged Gaussian integrations. The results vary the most at small 
frequency, and appear to systematically approach the exact re- 
sult as At — + 0. The results are less sensitive to At for larger 
frequencies, and since we have akeady seen that reducing the 
step size in e tends to only smooth out the oscillations (without 
changing the shape too much), the errors at higher frequencies 
must be coming from the finite cutoff tmax- 

Finally, we show the Green's function as a function of tid in 
Fig-El The imaginary part is in the main plot, and the real part is 
in the inset. In the exact result (blue curve), the imaginary part 
is an even function, and the real part is an odd function. The 
results of the nonequilibrium calculation do not share this sym- 
metry, but it appears to be getting restored as At 0. More 
problematic is the issue of the value of g"^ at trei = 0, which 
is determined solely by the electron filling. The results appear 
to be getting worse as At 0. It must be that if we increase 
the time cutoff, the results will ultimately start moving back to- 
ward their correct value, but we cannot check this explicitly due 
to the finite computer resources that are available. 

We can be more quantitative about the short time behavior 
though. To do so, we can calculate the first few moments of 
the integrals of the Green's function and the self-energy over cj, 
and compare them to the exact results for those integrals. The 



zeroth moment relates to the function at tioi = 0, the first power 
to the slope, and the second power to the curvature. When we 
examine the results for U — 0.5 and [/ = 1, we find the zeroth 
moment of g and E is in error by about 7%, the first moment 
by 10% for g and 20% for S, and the second moment of g 
and S by 15 — 20%. The results do not depend too strongly 
on the step size for e in the integrations, as expected, because 
the oscillations average out of the moments. The first moment 
appears to extrapolate to its correct value as At ^ for the 
A^ — 54, 55 Gaussian integration with a fixed tmax, but the 
zeroth and second moments do not appear as if they scale to 
the right result. This is most likely due to the fact that the t^ax 
needs to be increased, but the large step size for e may also play 
a role. 

When we try to examine the insulating phase, we find the 
agreement with the exact results becomes much worse. This is 
because there are low-energy features which require large times 
to be determined accurately. Also, the larger U , the smaller At 
needs to be to obtain good accuracy. Our results for J7 = 1.5 
are too preliminary to report quantitative values here. 

The important question is whether these low-energy features 
survive as the electric field is turned on. If they are destroyed 
by the field, then the computational scheme that we are using 
should be able to accurately determine nonequilibrium results 
at short times. If they survive, then it will be difficult to get 
high accuracy results for the nonequilibrium case in the Mott 
insulator. Since the presence of a field pumps energy into the 
system, and that energy can be used to create excitations across 
a gap, it is easy to believe that the gap features do not survive the 
introduction of a large electric field, but we cannot definitively 
say whether this is actually true at this point. 

V. Bloch oscillations 

One of the most interesting nonlinear phenomena of a mate- 
rial is the production of Bloch oscillations in the current as a 
function of time when a constant idc) electric field is applied. 
In the absence of interactions (which cause scattering), the cur- 
rent will oscillate forever, with a constant amplitude (the period 
is determined by the strength of the electric field). As we turn 
on the scattering, we expect the oscillations to be dampened, 
but perhaps to maintain the same period (and even survive in 
the steady state). If the scattering becomes large enough, then 
the oscillations should disappear completely. By calculating the 
Green's functions in the presence of a field that is turned on at 
r = 0, we can study how the current initially starts, and how 
it evolves into a steady state. Due to the need for a finite tmax, 
we can only go so far out in time before the calculation must 
terminate. 

In Fig. |6l we plot the Bloch oscillations of the current as a 
function of the average time, for U = 0.5, (3 = 1, E — 1, and 
the N = 54, 55 Gaussian integration scheme. The t,nax is equal 
to 15. As At is reduced, one can clearly see that the results are 
beginning to converge. Furthermore, it is also clear that there 
is a damping of the transient response as we move forward in 
time. We examine the behavior of this transient damping in 
Fig.El where we have preliminary results for the current density 
for three values of U. The results for U — 0.5 and [/ = 1, are 
calculated with a step size of At = 0.05. At the moment, it 



is not clear how much of an effect the boundaries at tniax have 
on the results, but it may be that the largest time results are not 
fully trustworthy. Some interesting behavior can be seen in the 
figure: (i) we see the damping that is expected, and that it gets 
more strongly damped as U is increased; (ii) as U is increased, 
it appears that the period may be decreasing, which is not an 
expected result; and (iii) the behavior of the transient evolution 
is quite complex, and we do not appear to have reached the 
steady state yet. 

VI. Conclusions 

In conclusion, we have shown that there is a straightforward 
way to perform many-body physics calculations in real time for 
both equilibrium and nonequilibrium situations by formulating 
the problem on a Kadanoff-Baym contour and discretizing it. 
There are a number of numerical issues that arise from this ap- 
proach, coming from the discretization of the contour and its 
truncation, as well as from the discretization of the energy space 
needed to perform numerical quadratures. By using the equihb- 
rium results as a benchmark, we can see how the different dis- 
cretization operations affect the overall accuracy of the calcula- 
tions. It seems like the small-time behavior can be understood 
fairly well by using reasonable choices for the discretizations 
and the time-domain cutoffs. It is more difficult to obtain good 
results for the longer-time behavior. One also needs to have 
good control of the numerics to be able to accurately perform 
Fourier transformations to real frequency. Maintaining control 
of these different approximations is made more difficult by the 
finite computational power that is available to solve these prob- 
lems. 

The fact that we must iterate our equations to a self- 
consistent solution brings a number of unknown issues to the 
table. First, we have modified the Green's function by intro- 
ducing a time-domain cutoff, which is artificially changing the 
boundary conditions. It is well known that Green's functions 
are determined uniquely by their equation of motion and their 
boundary conditions. How much of an effect the change of the 
boundary condition has on the results is difficult to estimate be- 
cause of the nonlinearities introduced by the iterative solution. 
Second, we are not able to iterate the equations for an infinitely 
long period of time, so the smaller we make the discretizations, 
the fewer iterations we are able to complete. For example, when 
the matrices have a size on the order of 700 x 700, we can eas- 
ily perform 50 or more iterations, but we are reduced to about 7 
iterations when they are 2000 x 2000. It is difficult to tell how 
much error is introduced by this. 

In the future, we will use this technique to study and analyze 
better the behavior of strongly correlated materials in a large 
electric field. Eventually, we hope to be able to generalize this 
approach to apply it to multilayered nanostructures and thereby 
be able to directly calculate the current- voltage characteristic of 
a strongly correlated device. 
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